{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculating rigorous bounds on $\\pi$ with interval arithmetic \n",
    "\n",
    "[David P. Sanders](http://sistemas.fciencias.unam.mx/~dsanders/)\n",
    "\n",
    "Department of Physics, Faculty of Sciences, National University of Mexico (UNAM)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following the book [*Validated Numerics*](http://press.princeton.edu/titles/9488.html) (Princeton, 2011) by [Warwick Tucker](http://www2.math.uu.se/~warwick/CAPA/warwick/warwick.html), we find *rigorous* (i.e., guaranteed, or *validated*) bounds on $\\pi$ using interval arithmetic, via the [`IntervalArithmetic.jl`](https://github.com/dpsanders/IntervalArithmetic.jl) package.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using IntervalArithmetic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating $\\pi$ via a simple sum "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many ways to calculate $\\pi$. For illustrative purposes, we will use the following sum\n",
    "\n",
    "$$ S := \\sum_{n=1}^\\infty \\frac{1}{n^2}.$$\n",
    "\n",
    "It is [known](http://en.wikipedia.org/wiki/Basel_problem) that the exact value is $S = \\frac{\\pi^2}{6}$.\n",
    "Thus, if we can calculate rigorous bounds on $S$, then we can find rigorous bounds on $\\pi$. \n",
    "\n",
    "The idea is to split $S$ up into two parts, $S = S_N + T_N$, with\n",
    "$ S_N := \\sum_{n=1}^N \\frac{1}{n^2}$\n",
    "and $T_N := S - S_N = \\sum_{n=N+1}^\\infty n^{-2}$.\n",
    "\n",
    "By bounding $T_N$ using integrals from below and above, we can see that $\\frac{1}{N+1} \\le T_N \\le \\frac{1}{N}$.\n",
    "Rigorous bounds on $S_N$ are found by calculating it using interval arithmetic.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$S_N$ may be calculated by summing either forwards or backwards. A naive (non-interval) version could be the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "forward_sum_naive (generic function with 1 method)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function forward_sum_naive(N)\n",
    "    S_N = 0.0\n",
    "\n",
    "    for i in 1:N\n",
    "        S_N += 1./(i^2)\n",
    "    end\n",
    "\n",
    "    S_N\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.6448340718480652,9.999500016122376e-5)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S = forward_sum_naive(10000)\n",
    "err = abs(S - pi^2/6.)  # error\n",
    "S, err"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interval arithmetic "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find *rigorous* bounds for $S_N$, we use interval arithmetic: each term is enclosed in an interval that is guaranteed to contain the true real value. A first idea is simply to wrap each term in the `@interval` macro, which converts its arguments into containing intervals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "forward_S_N (generic function with 1 method)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function forward_S_N(N)\n",
    "    S_N = @interval(0.0)\n",
    "\n",
    "    for i in 1:N\n",
    "        S_N += @interval( 1./(i^2) )\n",
    "    end\n",
    "\n",
    "    S_N\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "elapsed time: 0.787551828 seconds (236886656 bytes allocated, 24.83% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[1.6449240668871583, 1.644924066909359]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 10^5\n",
    "@time rigorous_approx_S_N = forward_S_N(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We incorporate the respective bound on $T_N$ to obtain the bounds on $S$, and hence on $\\pi$. We can also optimize the code by creating the interval `@interval(1)` only once:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "elapsed time: 4.652939682 seconds (2080127080 bytes allocated, 27.31% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "([3.1415926534833463, 3.1415926536963346],2.1298829366855898e-10)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function forward_sum(N)\n",
    "    S_N = @interval(0.0)\n",
    "    interval_one = @interval(1.)\n",
    "\n",
    "    for i in 1:N\n",
    "        S_N += interval_one / (i^2)\n",
    "    end\n",
    "    \n",
    "    T_N = interval_one / @interval(N, N+1)\n",
    "    S = S_N + T_N\n",
    "\n",
    "    sqrt(6*S)\n",
    "end\n",
    "\n",
    "N = 10^6\n",
    "@time S = forward_sum(N)\n",
    "S, diam(S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can ask for the midpoint--radius representation, which shows that the calculated bounds are correct to around 10 decimal places:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3.141592653589824,5.834666083615048e-11)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "midpoint_radius(S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We may check that the true value of $\\pi$ is indeed contained in the interval:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "big(pi) ∈ S"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We may repeat the calculation, but now summing in the opposite direction. Due to the way that floating-point arithmetic works, this gives a more accurate answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "elapsed time: 5.638430487 seconds (2080153808 bytes allocated, 28.92% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "([3.1415926535893144, 3.141592653590272],9.57456336436735e-13)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function reverse_sum(N)\n",
    "    S_N = @interval(0.0)\n",
    "    interval_one = @interval(1.)\n",
    "\n",
    "    for i in N:-1:1\n",
    "        S_N += interval_one / (i^2)\n",
    "    end\n",
    "    \n",
    "    T_N = interval_one / @interval(N, N+1)\n",
    "    S = S_N + T_N\n",
    "\n",
    "    sqrt(6*S)\n",
    "end\n",
    "\n",
    "N = 10^6\n",
    "@time S = reverse_sum(N)\n",
    "S, diam(S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the `sqrt` function is guaranteed (by the [IEEE 754 standard](http://en.wikipedia.org/wiki/IEEE_floating_point)) to give correctly-rounded results, so the resulting bounds on $\\pi$ are *guaranteed to be correct*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note also that due to the way that the `IntervalArithmetic` package works, we can make the code simpler, at the expense of some performance. Only two explicit calls to the `@interval` macro are now required:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "elapsed time: 6.699738655 seconds (2264155172 bytes allocated, 27.45% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "([3.1415926535893144, 3.141592653590272],9.57456336436735e-13)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function reverse_sum2(N)\n",
    "    S_N = @interval(0.0)\n",
    "\n",
    "    for i in N:-1:1\n",
    "        S_N += 1 / (i^2)\n",
    "    end\n",
    "    \n",
    "    T_N = 1 / @interval(N, N+1)\n",
    "    S = S_N + T_N\n",
    "\n",
    "    sqrt(6*S)\n",
    "end\n",
    "\n",
    "N = 10^6\n",
    "@time S = reverse_sum2(N)\n",
    "S, diam(S)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.3.8-pre",
   "language": "julia",
   "name": "julia 0.3"
  },
  "language_info": {
   "name": "julia",
   "version": "0.3.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
